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Abstract We present the first method to handle cur- 
vature regularity in region-based image segmentation 
and inpainting that is independent of initialization. 

To this end we start from a new formulation of 
length-based optimization schemes, based on surface 
continuation constraints, and discuss the connections 
to existing schemes. The formulation is based on a cell 
complex and considers basic regions and boundary ele- 
ments. The corresponding optimization problem is cast 
as an integer linear program. 

We then show how the method can be extended to 
include curvature regularity, again cast as an integer lin- 
ear program. Here, we are considering pairs of boundary 
elements to reflect curvature. Moreover, a constraint set 
is derived to ensure that the boundary variables indeed 
reflect the boundary of the regions described by the 
region variables. 

We show that by solving the linear programming re- 
laxation one gets quite close to the global optimum, and 
that curvature regularity is indeed much better suited 
in the presence of long and thin objects compared to 
standard length regularity. 

1 Introduction 

Regularization is of central importance for many in- 
verse problems in computer vision including image seg- 
mentation and inpainting [19,9, 18,2,29,4 . The intro- 
duction of higher-order regularizers in respective en- 
ergy minimization approaches is known to give rise to 
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substantial computational challenges. Some of the most 
powerful approaches to image segmentation are based 
on region integrals with regularity terms defined on 
the region boundaries [3| fT9| f2T | f5|t9|fT2 | fT6] . While many 
such methods make use of length as a regularity term, 
only few use curvature regularity. This is in contrast 
to psychophysical experiments on contour completion 
[T5] where curvature was identified as a vital part of 
human perception. Note that curvature regularity is 
qualitatively different from length regularity. As shorter 
boundary curves are preferred in the length case, this 
causes the well-known shrinking bias. This is not the 
case for curvature, since the total curvature of any 
closed, convex curve is equal to 2ir. 

Length regularization has become an established 
paradigm as there exist many powerful algorithms for 
computing optimal solutions for such energy function- 
al, either using discrete graph-theoretic approaches 
based on the min-cut /max- flow duality [14,5 or us- 
ing continuous PDE-based approaches using convex re- 
laxation and thresholding theorems [20]. Region-based 
problems for segmentation using curvature regularity 
have typically been optimized using local optimization 
methods only (cf. [2T | [T2 j ). As a consequence, experi- 
mental results highly depend on the choice of initializa- 
tion. Moreover, these methods do not offer any insights 
concerning how close the computed solution is to the 
(unknown) global solution. 

In this paper, we propose a relaxed version of region- 
based segmentation which can be solved optimally. The 
key idea is to cast the problem of region-segmentation 
with curvature regularity as an integer linear program 
(ILP). By solving its LP-relaxation and thresholding 
the solution we obtain a solution to the original integer 
problem and are able to evaluate a bound on its quality 
with respect to the globally optimal solution. In addi- 
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thresholding scheme 
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Fig. 1 With the proposed method, long and thin structures are much better handled than with length-based approaches. 




damaged image with length regularity 

Fig. 2 Curvature regularity on the level lines improves inpainting. 



with curvature regularity 



tion, we show that the method readily extends to the 
problem of inpainting. 

Figure [T] demonstrates that the proposed method 
allows segmenting objects in a way which preserves 
perceptually important thin and elongated parts. The 
found solution is within 1.3% of the global optimum. 
Figure [2] demonstrates the superior performance of cur- 
vature regularity over length regularity in a correspond- 
ing inpainting experiment. 

Existing Work on Curvature Regularity. For contour- 
or edge-based segmentation methods researchers have 
successfully developed algorithms to optimally impose 
curvature regularity using shortest path approaches pQ 
or ratio cycle formulations [24] on a graph representing 
the product space of image pixels and tangent angles 
[22J. In the region-based settings considered, curvature 
is usually handled by local evolution methods [8,12,21, 
[29] . Among the methods that pre-date our conference 
publication [25], the only exception we are aware of is 
the inpainting approach of Masnou and Morel [18] who 
can optimize the Li-norm of the curvature in the ab- 
sence of regional data terms using dynamic program- 
ming. 

In this paper we propose an LP-relaxation approach 
to minimize curvature in region-based settings. In con- 
trast to [18 it allows imposing arbitrary functions of 



curvature and arbitrary data terms. The algorithmic 
formulation is based on the concepts of cell complexes 
and surface continuation constraints which have been 
pioneered by Sullivan [28^ and Grady [13] in the context 
of 3D-surface completion. 

The present paper is based on our preliminary work 
in [25], but contains several novelties. Firstly, we show 
that the constraint system in [25] needs to be aug- 
mented by additional constraints in order to ensure that 
the boundary of the region-based segmentation is cor- 
rectly estimated. Secondly, we discuss the connections 
of our method, when restricted to length regularity, 
with standard methods for length regularity and com- 
pare experimentally. In addition, the original inpainting 
method has improved by incorporating a boundary es- 
timation scheme. 

There has been a subsequent paper by El-Zehiry and 
Grady [11] which optimizes the same model as in our 
original paper, but applies quadratic pseudo-Boolean 
optimization (QPBO) for obtaining the solution. In the 
case of a square grid, QPBO is able to efficiently com- 
pute a solution with curvature regularization, and with 
a subsequent probing stage often the global optimum is 
found. However, the discretization artefacts are severe 
for a cell complex of squares. Better connectivities can 
only be handled approximately. 
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The source code associated to this paper is freely 
available at 
http : //www. maths . 1th. se/matematiklth/personal/ 
tosch/ download . html. 



2 Length-based Segmentation Problems 

To detail the proposed method for curvature regular- 
ity, we first introduce a novel method for length-based 
segmentation problems. In practice there are more effi- 
cient algorithms for this problem [6 ,20 , but in contrast 
to them the presented one is easily extended to curva- 
ture. A comparison to the known techniques is given at 
the end of this section. 

Given an image I : 4? — >• the problem is to seg- 
ment it into two regions, foreground and background. 
Here "region" means an arbitrary subset of i?, i.e. there 
can be several disconnected components and each one 
can have holes. Hence, each point x G Q is to be as- 
signed a region iz(x) G {0, 1} where denotes back- 
ground, 1 foreground. 

The desired segmentation is defined as the global 
optimum of an energy function, consisting of two terms. 
The first one is called the data term and specified by a 
function #o( x ) for points belonging to the background 
and a function #i(x) for the foreground. Both functions 
will generally depend on the input image I. In addition 
there is a regularity term that penalizes the length of 
the segmentation boundary by a weighting parameter 
v > 0. The arising energy minimization problem to be 
solved is 
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-u\C\ 



(1) 



where C = <9{x | Viz(x) = 0} is the set of points where 
u is discontinuous (i.e. "jumps" from to 1) and \C\ 
denotes its one-dimensional measure. In other words, 
C is a set of closed lines (which may include parts of 
the boundary of Q) and \C\ denotes the sum of the 
length of all lines. 

For convenience, we reformulate ([I]) by splitting the 
integrand of the first term into a constant term and 
one depending on u. Defining g{x) = #i(x) — #o( x ) the 
resulting functional is 



mm 



g(x.) u(x.) -\-v\C\ + const . 



(2) 



In the following we will ignore the constant except when 
evaluating relative gaps between a lower bound and the 
energy of some segmentation. 




Fig. 3 The basic concepts of our method, (a) A cell complex, 
(b) The method considers oriented regions and oriented boundary 
elements. 



2.1 Discretization 

In this paper we consider discretized segmentation 
problems where instead of optimizing infinitely many 
values {u(x) | x G 42} we only consider finitely many 
"basic regions" and jointly assign all points in a basic 
region to the same segment. Note that in practice we 
will always get a discrete input image /, where the ba- 
sic regions are given by pixels. Hence, the data term 
by itself will produce such an assignment even for the 
continuous problem. This is no longer true when the 
regularity term is added, but in practice the discretized 
energy function can be designed to account for this phe- 
nomenon 16II20]. 

We require that our set of basic regions - denoted 
T - be a cell complex and a partitioning of 42, i.e. that 
(1) no two regions overlap and (2) the union of all basic 
regions yields 42. An example is given in Figure [3] (a). 

The presented approach makes use of another essen- 
tial part of a cell complex: boundary segments. These 
are the line segments that form the borders of the basic 
regions. Usually a boundary segment has two neigh- 
boring regions, except for segments at the border of 42 
where there is only one. The set of all boundary seg- 
ments is denoted £ . As will be shown below we need to 
consider both possible ways of traversing a boundary 
segment. Hence, for each boundary segment we con- 
sider two oriented boundary segments, also called line 
segments in the following. The set of all these segments 
is denoted £° and t{e) will denote the length of a line 
segment e. 

In essence, the data term will be defined in terms of 
the basic regions, the regularity term in terms of bound- 
ary segments. To approximate the continuous problem 
sufficiently well, basic regions should generally not cor- 
respond to pixels. Instead, as detailed in Figure [4] we 
split the pixels into either 4 or 32 basic regions, in- 
duced by lines with 8 and 16 different directions re- 
spectively. In the following we will refer to this as 8- 
and 16-connectivity. 
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^-connectivity 



16-connectivity 



Fig. 4 Splitting a pixel into basic regions using lines with 8 and 
16 different directions. 



2.2 An Integer Linear Program 

The presented method casts a discretized version of ([!]) 
as a so-called integer linear program, i.e. minimizing a 
linear cost function over integral variables and subject 
to linear constraints. There are two sets of variables: 
firstly, for each basic region / G T there is a region 
variable G {0, 1} with indicating that the region 
belongs to the background, 1 to the foreground. The 
second set contains a boundary variable y B G {0, 1} 
for every oriented boundary segment e G £° . Here we 
want yf to be 1 only if exactly one of the adjacent basic 
regions belongs to the foreground. Now we are already 
in a position to express the cost function: 



(3) 



where c# contains entries 



g(x) dx 



Formalizing these constraints (one for each bound- 
ary segment) involves the concept of orientations for 
both regions and boundary segments. For boundaries 
we have already introduced this concept, but it is es- 
sential for the constraint system that we (arbitrarily) 
define a "positive" and a "negative" orientation for each 
boundary segment. 

For a region, an orientation denotes one of the two 
possibilities for traversing its boundary line - clockwise 
or counter-clockwise. Here it is essential that all regions 
have the same orientation. 

Now, to formalize the surface continuation con- 
straint we define the notion of positive and negative 
incidence of regions / G T and line segments I G S° 
to boundary segments e G £ . Ultimately the constraint 
will then state that the weighted sum of "active" re- 
gion incidences must be equal to the weighted sum of 
"active" boundary incidences, where active refers to el- 
ements where the associated indicator variable is 1. 

The incidence of a region / to a boundary segment 
e is denoted m{ (for "match") and defined as unless 
/ contains e in its boundary line. Otherwise, m{ is 1 
if according to the orientation of / the segment e is 
traversed in its positive orientation, otherwise —1 for 
the negative orientation. The incidence of a line I and 
a boundary segment e is denoted m l e and defined as 
unless I is an orientation of e, otherwise 1 if I is the 
positive and — 1 if I is the negative orientation of e. The 
constraint system can now be stated as 
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For clarity, let us look at the following example: 



and Cb contains entries c e B = vt(e). Note that since 
in our case basic regions are always subsets of a single 
pixel this weight is simply the function g(-) evaluated 
at the pixel times the area of the basic region. 

Since v is positive minimizing ([3| by itself would 
set all boundary variables to 0, so we need constraints. 
Indeed, up to a few ambiguities (see next section) the 
region boundary is completely specified by the region 
variables and the boundary variables serve to render 
the cost function linear. They are forced to describe 
the correct boundary by the linear constraint system 
that we now describe. In words it can be stated as 

Surface Continuation Constraint: When- 
ever a basic region is part of the foreground, 
along each of its boundary segments the fore- 
ground must either continue with another fore- 
ground region or with an appropriately oriented 
boundary segment. 




Here the constraint for the bold edge reads 



y A R -v B R 



(5) 



In the case where = 1 and yj^ = 0, i.e. region A 
is foreground and region B is background, y l B will be 
forced to 1, whereas y l B will be set to 0. If instead B is 
foreground and A background, this will force y l B to be 
1 and y l B to be 0. If both A and B belong to the same 
component, the constraint leaves some freedom for the 
boundary variables: they can now both be or both 
be 1. The latter is undesirable, but will not happen as 
long as the length weight is strictly positive. However, 
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when we integrate curvature below we will need extra 
constraints to prevent this case. 

Finally, we summarize the integer linear program to 
be solved as: 



min c^y^ + y B 
y 



(6) 



y£e{0,i} V/G.F 

y l B G{0,1} Vl€£° . 

As we will show now this problem can be efficiently 
solved by computing a graph min-cut. 




Standard Graph Cut 



Subdivision 



Fig. 5 Comparison of standard graph cuts and the novel sub- 
division scheme, both run with the same length weight and an 
8-connectivity. (Input image shown in Figure [TT] ) 



2.3 Relation to Graph Cuts 

Discrete approaches to image segmentation are very 
well studied in computer vision and the vast major- 
ity uses pixels as their basic regions. However, they do 
usually not express (length) regularity in terms of the 
boundary segments of these regions - this would imply 
a four- connectivity. 

Instead, these approaches are based on graphs where 
the basic regions correspond to nodes and the length 
term is represented in terms of edges that connect pairs 
of nodes [6]. We denote the set of nodes V, where each 
p G V represents the center point of a basic region 
/ G T . Analogous to the above integer program, each 
center p of a basic region is associated a binary variable 
y p G {0, 1} indicating foreground and background. The 
smoothness term is modeled by a set of edges M (also 
called neighborhood) and for length regularity can be 
expressed as: 

lcl *m ( ^¥^\\ {1 - SM) ' 

{p,q)EM 

where k(M) is a normalization constant that ensures 
that the weights remain comparable when the size of 
the neighborhood is enlarged. 

Minimizing this energy can be written as a mini- 
mum cut problem, which again can be written as 



ypG{o,i}i^i 



v P ,q \y v -y q \ 



It is well-known, e.g. [10], that such absolutes can be 
rewritten as linear programs, i.e. this problem can be 
equivalently written as 

y min c T R y P + ^ w p , q (a+ g + a~ q ) 



S.t. X f) X ri 



p,q p,q 



%m V(p, q)eJ\f 



If we assume that the neighborhood links exactly 
all pairs of neighboring basic regions (i.e. those that 
share a boundary segment), these are exactly our sur- 
face continuation constraints ([5|. Since graph cuts can 
be optimized globally efficiently it follows that ^ is 
polynomial-time solvable. 

In summary, what we have proposed so far is a re- 
striction of graph cuts to graphs that are planar when 
source and sink are removed. This restriction is crucial 
for (our solution of) the problem we really want to solve: 
to integrate curvature regularity into the framework. 

In contrast to standard applications of graph cuts 
our method relies on subdiving pixels. In the case of a 
very strong data term a boundary pixel will probably 
not be split into a foreground and a background part. 
Standard graph cuts might then be the better way to 
reflect the length regularity term. However, in practice 
boundary pixels usually have weak data terms due to 
partial overlap. Figure [5] shows the resulting segmenta- 
tions of both schemes, where the input image can be 



y P G {0, , a+ q , a~ q > V(p, q) G M . 



found in Figure 11 Here it can be seen that the novel 
scheme gives access to a higher resolution and produces 
slightly different segmentations. These are actually of 
lower energy than for the standard scheme. Hence, if 
anything we have gained something with the novel for- 
mulation. 



2.4 Relation to Sullivan's Method 

Sullivan j27] l28] gave a method for finding a surface of 
minimal (weighted) area in 3D-spac^] subject to the 
constraint that the surface span a given (set of) bound- 
ary lines. His method also relies on a cell complex and 
can be written as a polynomial-time solvable minimum 
cost flow problem. Further, the method allows integrat- 
ing volume terms. 

1 In fact, he considered general N-D spaces, but this is not 
important here. 
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In essence, what we have done so far is a restriction 
of this method to 2D and where the set of prescribed 
boundary elements is empty (these objects would be 
one-dimensional). However, there is one major differ- 
ence: in Sullivan's method all variables are unrestricted, 
whereas we have the constraints y^ R G {0, 1} and y l B G 
{0, 1}. And of course our main goal is to integrate cur- 
vature regularity. 

3 Handling Curvature Regularity 

We now show how the described integer linear program 
can be generalized to curvature regularity, i.e. to the 
model 

min / Gf(x) u(x) + v \C\ 

u:fl->{0,l} J 
Q 

+\j Mx)|*W(x) . (7) 

c 

Here A > is a curvature weight, ftc( x ) stands for the 
curvature of the line C at a given point x of the line 
and p > is an arbitrary exponent (usually p = 2 is 
a good choice). The notation dH 1 (x) signifies that the 
integral is over a set of lines and that it is independent 
of the parameterization of these lines. 

We emphasize that our method allows more general 
regularity terms, namely arbitrary positive functions 
depending on position, direction and absolute curva- 
ture. In particular this allows spatially weighted regu- 
larity terms. 

3.1 Discretizing the Problem 

Again, our solution is based on a cell-complex and the 
data term is handled in the exact same way as above. 
That is, we again have region indicator variables Xf for 
all / G T . We were able to express the length regu- 
larity in terms of (single) boundary segments. For cur- 
vature, this is not possible: all boundary segments are 
straight lines, hence have curvature everywhere. The 
only points where non-zero curvature can occur are the 
meeting points of two boundary segments. Hence it is 
common to consider pairs of line segments [22| fTj to ex- 
press curvature regularity. So far, however, this was not 
compatible with region terms. 

For every pair of adjacent line segments with 
compatible orientations we now have an indicator vari- 
able y l g l2 G {0, 1}. Here we follow the convention that 
the line segment that is traversed earlier is also listed 
first in the pair. We get a cost function of the form 

c rYr + c bYb , (8) 




Fig. 6 The angle 6 is the basis for computing the curvature of 
the pair of black lines. 

where y# now contains all the pairwise variables. We 
proceed to describe the entries of the corresponding cost 
vector c|J. 

3.2 Computing the Weights 

Computing curvature from two adjacent line segments 
is based on considering the direction change, measured 
by the angle in Figure [6j 

There are basically two ways to compute the term 
\k\ p from this angle: firstly, we can just take the power p 
of the angle (which should be measured in arc length). 
The second method is based on the work of Bruckstein 
et al. [7 , and makes also use of the lengths and 
£(l2) of the two lines: 

^^' U»iw) )' ■ 

The original idea of Bruckstein et al. was to take the 
longest straight lines pre- and succeeding a direction 
change. In our context we can only consider elemen- 
tary boundary segments, so we do not have the same 
convergence properties. In practice we found that both 
weights work fine and Bruckstein et al's weights signif- 
icantly reduce the running time of the employed linear 
programming solver. 

Denoting either of the two arising weights as wi ly i 2 , 
we get one of two components for the cost entry Q l5 / 2 . 
Together with length regularity this entry is 

There are however some special cases involving the im- 
age border, a point that we neglected in [25]: firstly, at 
the four corners of a rectangular domain Q we will want 
to set the curvature weight to 0. Note that we still pe- 
nalize the angle in which a region boundary meets the 
domain boundary dfl. 

A second point relates to the length weight: when- 
ever one (or two) of the line segments in a pair is part of 
the domain border, its length should be set to 0. Oth- 
erwise there would be a bias towards associating the 
regions at the image border to the background. 
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3.3 An Adequate Constraint System 

As before for length regularity, the linear cost func- 
tion ^ needs to be minimized subject to suitable con- 
straints to closely reflect the discretized model function. 
In [25] we presented two sets of constraints, with the 
aim to ensure that the boundary variables indeed de- 
scribe a boundary of the region variables. In this work 
we show that two more sets of constraints are needed, 
where the second one becomes necessary only when sev- 
eral regions meet in a single point. In particular, we will 
show that in this case there are several valid boundaries 
and that our method searches for the one with least 
cost. 

Firstly, we adapt the surface continuation con- 
straints to the new kind of boundary variables. To this 
end, we define the incidence m l e ljl2 of a line segment 
pair and a boundary segment e G S . This value is de- 
fined as unless l\ is an orientation of e (note that I2 
is irrelevant for this value). Otherwise the value is the 
same as the previously defined incidence of l\ and e. 
Hence, the surface continuation constraints read 



f f 
mi y J R 



£ 



Ve e S 



(9) 



These constraints alone leave a lot of freedom. In par- 
ticular, we could choose line pairs without direction 
changes everywhere. What is really wanted, though, is 
that for every active y l b 1,Ii2 both l\ and I2 belong to the 
region boundary induced by the region variables. This 
is ensured by two sets of constraints, where the first is 
called boundary continuation. In words it can be stated 
as 

Boundary Continuation Constraint: If a 

pair of line segments h,h is active, there must 
be a succeeding pairl2, 13 that is also active. Like- 
wise there must be a preceeding active pair /o^i- 

These constraints ensure that the active line pairs actu- 
ally define closed paths. They are identical to the con- 
straints arising for the computation of shortest paths 
in a graph such as [I] and are stated as 

E^=E^ 2 Vhe£° . 

Now we have paths, but we cannot guarantee that all 
parts of these paths are actually region boundaries. In- 
deed, we invite the reader to check that the configura- 
tion in Figure [7] (a) satisfies all constraints introduced 
so far. Moreover, for small length weights and squared 
curvature the cost of this configuration will be lower 
than those of the desired configuration shown in part 
(b). To exclude cases such as (a) from the optimization, 
we add a new set of constraints: 






(a) 



(b) 



Fig. 7 Without the boundary consistency constraints the con- 
figuration in (a) is valid and - for squared curvature and with- 
out length penalty- cheaper than the desired one in (b). With 
boundary consistency (b) remains feasible, but (a) is excluded, 
as desired. 

\1I 



(a) 



(b) 



(c) 



Fig. 8 Determining boundaries is ambiguous: (a) A segmenta- 
tion where two regions meet in a point, (b) A low-cost boundary 
configuration, where different pairs are indicated by different col- 
ors, (c) An equally valid, but more expensive boundary. 



Boundary Consistency Constraint: For ev- 
ery boundary segment, only one of the two pos- 
sible orientations can be active. 

To formalize this, we denote by e~^ and e^~ the positive 
and negative orientation of a boundary segment e. Then 
the constraint can be written as 

E^ + E^f' 2 ^ 1 Ve ^ • 

Similar sets of constraints can be derived, e.g. 
Ez x Vb A + T,i 2 V l B e < but experimentally we 
found them to be redundant. Moreover, if e is part of 
the domain border df2 then one of the orientations will 
never occur in a desired configuration. We simply dis- 
regard the corresponding variables. 

With these constraints, there is still one issue left 
to take care of, and it affects cases where several re- 
gions meet in a point. In this case, as shown in Figure 
[5J there are several valid boundary configurations. As 
long as only two regions meet, the given constraint sys- 
tem ensures that indeed the configuration with lower 
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Fig. 9 If three or more regions meet in a point, self-intersecting 
boundaries can define valid segmentations. If these are undesired, 
a fourth constraint set is needed. The different colors denote dif- 
ferent phases when traversing the line. There is only one line. 



cost is selected. However, as soon as there are three or 
more regions meeting in a point, this constraint system 
will allow configurations with crossings, as exemplified 
in Figure [9] Should we really avoid such configurations? 
There is actually no obvious answer if one has in mind 
that, in the theory of continuous plane curves, the set 
enclosed by a positively oriented closed curve can be de- 
fined as the points of positive index (also called winding 
number) [23 . In the example of Figure [9j the index of 
each point in a black triangle with respect to the outer 
curve defined by the arrows is one, so it makes sense 
to consider the curve as a region boundary For the 
sake of comparison and completeness, we provide in sec- 
tion [6| Figure [13] a couple of experiments done either 
with or without crossing prevention. In the former case, 
we add a new set of constraints: 



Crossing Prevention Constraint: If two 

pairs of line segments cross, only one of them 
may be active. 

Denoting C the set of crossing line pairs, this con- 
straint is easily formalized as 



Vb 



y l i M <i V(Zi,Z 2 ,Z 3 ,Z 4 )eC . 



In practice these constraints are ignored in a first phase. 
Afterwards, the (usually very few) violated constraints 
are added in passes and the system is re-solved, until 
there are no more violated constraints. In our experi- 
ments we never needed more than 9 passes. However, 
solving the arising programs can be quite time consum- 
ing, even when starting from the previous configuration. 

In summary, region-based segmentation with curva- 
ture regularity is now expressed as the integer linear 

2 the index of a point on a discrete grid with re- 
spect to a discrete curve can also be computed, 
see for instance the winding number algorithm at 
www.softsurfer.com/Archive/algorithm_0103/algorithm_0103.htm 



program 

min c^y R + y B 

YR,yB 

S.t. £>effl = m e' l2 yB h Vee£ 

f€T h,he£° 

!/^{0,l}V/eJ, y%' 12 € {0,1} V h,l 2 G £° , 
and the optional constraints 

y h,h +y h,u< 1 v(Z 1 ,Z 2 ,Z 3 ,y eC. 

A first evaluation of this new scheme was given in 
Figure [l] on page [2] which shows that curvature regular- 
ity is better suited to ensure connected regions in the 
presence of long and thin objects. Note that the found 
solutions for curvature are generally not globally opti- 
mal. Details on the optimization scheme are given in 
Section [5] and more experiments will be given in Sec- 
tion [6] First, however, we discuss how to handle the 
problem of image inpainting. 

4 Inpainting 

In image inpainting we are given an image I : f2 M 
together with a damaged region C Q. This damaged 
region can have arbitrarily many connected components 
and each of these can enclose holes. The task is to fill 
the damaged region with values that fit nicely with the 
values of / outside the damaged region. 

To this end, the above integer linear program is gen- 
eralized to give a structured inpainting approach, where 
we consider the continuous model [T8 l[T7P 8] 



,:Q- 



nn i} | J \K Cu ^)\ P dU\*)dt 



(10) 



s.t. u(x) = j(x) Vx e n \ n d 



where I\ and I u are the minimal and maximal intensities 
of / along the border of the damaged region, C Ujt = 
{x I u(x) = t} is the set of level lines for level t of u(-) 
and again p > is an exponent for curvature. 

For the case of absolute curvature (p = 1) and that 
Qd consists of 5m^/^-connected components (i.e. the 
components do not enclose holes), an efficient global 
optimization scheme was given in [18 , 17 . We are inter- 
ested in the more general problem of arbitrary domains 
and exponents p > 0. In particular, p = 2 is usually a 
better model. 
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4.1 Discretization 



As before for image segmentation, our strategy is to 



discretize the model (10) by introducing basic regions 
and pairs of boundary segments. However, the regions 
are no longer limited to two labels: the intensity Uf of 
region / G T can be anywhere between // and I u . We 
follow the strategy in [T8| fT7] and consider only integral 
values inside this range. The result is a fully discrete 
labeling problem. 

Naturally, one also needs to change the right-hand- 
side values of the inequality constraints. Moreover, for 
inpainting we always include the crossing prevention 
constraints since by definition level lines cannot cross. 




Fig. 10 Best viewed in color. Estimated boundary variables for 
the given intensity profiles. Colors denote different pairs (not all 
pairs are shown) (a) With the proposed approximation we do 
not always get level lines. Here the boundary variables may jump 
between different intensity levels, (b) With level variables one 
would get the correct solution, but the computational demands 
are too high in practice. 



To make sure that the boundary variables truly 
reflect level lines, it would be advisable to associate 
each basic region multiple variables, reflecting level sets. 
Then, there would be a binary variable for every in- 
tegral value I\ < k < I u where a value of 1 reflects that 
the intensity Uf of the basic region is at least equal to 
k (i.e. Uf > k). For k' > k this would naturally en- 
tail the constraints y k > y k . Moreover, there would 
be binary variables y* j that would be forced to be 
consistent with the level variables exactly as for binary 
segmentation. 

This strategy is however not practicable for the do- 
main sizes we want to address: the problem is way 
too large scale. Instead, we settle for one integral vari- 
able y f R G {//,-.., Iu] for every region / G di- 
rectly reflecting the intensity Uf of the face. In ad- 
dition, there are boundary variables reflecting the in- 
tensity differences between neighboring regions. Inside 
the damaged domain they can be restricted to values 
y l g l2 £ [0? I u ~ h\- For the fixed part, we only consider 
basic regions that border on the damaged region. At 
their other borders the respective boundary variables 
can take values in {0, . . . , I u }. In practice it is advis- 
able to first subtract the constant I\ from the entire 
image (note that each connected component of can 
be processed independently). 



The arising integer linear program is stated as 



As showrj^] in Figure 10 this strategy does generally 
not give level lines, but it is a reasonable approximation. 



mm c B y B 

Yr,Yb 



(ii) 



S.t. Ys^eVR = E m e ,h yB h Ve € £ 



o 



h h 

y h,h +y w*<j u v(i lt i 2 ,i 3 ,i A )eC 

y f R = I f -I l Vfcf2\f2 d 

y£e{0,... ,!„-/,} V/Cf2 d 

y% M e {0,..., /„-/,} Vh,l 2 e£° . 

Note that we have one such program for every con- 
nected component of 



4.2 Estimating Incoming Level Lines 

A correct estimation of the direction of the level lines 
touching from outside the inpainting domain i? is im- 
portant, since the method aims at prolonging them with 
no additional curvature, if possible. We follow the sim- 
ple and efficient method proposed by Bornemann and 
Maxz in [4 after the work of Weickert [30,31 on the 
robust determination of coherence directions in image: 
at a point x G Q\Q^ the coherence direction is the 
normalized eigenvector associated to the minimal eigen- 
value of the structure tensor @j: 



Many thanks to Yubin Kuang for providing these images. 



J(x) 



ln\n d V7 ff ®V7, 



(k p * l n \n d )(x) 



(12) 
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where ^Q\Q d is the characteristic function of f2\f2d, I a 
is defined as 



Krr * 



K n *t 



0\Q d 



(13) 



and K pi K a are Gaussian smoothing kernels with stan- 
dard deviations p and a. Experimentally, setting <j = 
1.5, p = 4 yields a reliable estimation of the incoming 
level lines directions along dfid- 



5 Optimization Strategies 

In general, solving integer linear programs is an NP- 
hard problem [26]. In some cases, where there are lin- 
ear inequality constraints Ax < b and the matrix A 
is totally unimodular one can find the global optimum 
by solving the linear programming relaxation [26 , i.e. 
the problem one obtains when dropping all integrality 
constraints on the variables. For instance, a constraint 
Hi £ {0,1} will be relaxed to yi £ [0,1]. The arising 
problem can be solved in (weakly) polynomial time us- 
ing interior point methods [32] . 

Several of the discussed systems are in fact 
polynomial-time solvable, in particular the length- 
based problems and the boundary continuation con- 
straints by themselves. The integer linear program for 
curvature regularity is however not in this class. Still, 
experimentally we found that solving the linear pro- 
gramming relaxation often gives nearly integral solu- 
tions. Moreover, the relaxation value provides a (usu- 
ally quite tight) lower bound on the original problem. 
We proceed to discuss this strategy in detail. 



5.1 Solving the Linear Programming Relaxation 

There are two popular ways of solving general linear 
programs. Firstly, there is the dual simplex method 
[TO] , based on refactorizing the constraint system. It is 
usually the most memory saving method and in prac- 
tice superior to the primal simplex method. For some 
variants of the simplex method exponential worst-case 
run-times have been proved. On practical problems the 
method often works very well, and there are competi- 
tive and freely available implementations, e.g. the solver 
Clp0 We found it quite useful for an 8-connectivity, but 
for a 16-connectivity - and very low length weights - we 
got acceptable running times only for some of the im- 
ages we tried. In other cases the solver was terminated 
after several days without having solved the problem. 
In the end we chose the length weights high enough to 



get acceptable running times. The results we get indi- 
cate that these settings actually provide a very good 
model. 

On the other hand, there are interior point methods 
which usually perform Newton-iterations on a primal- 
dual formulation of the problem. This entails frequent 
matrix inversions, solved via the sparse Cholesky de- 
composition. We are not aware of a freely available 
solver that performs well on large scale problems such 
as ours. Hence, we tested several commercial packages 
and found Gurobi and FICO Xpress to be well-suited 
solvers for our problem. Generally the running times of 
commercial interior point solvers are quite predictable. 
The downside is the memory consumption: we found 
that a little more than twice the memory consumption 
of a simplex solver is needed for our problem, where we 
tested with an 8- and a 16-connectivity. 

For segmentation, the combination of licensing is- 
sues and the high memory demands made us use the 
simplex method. For inpainting, where the memory de- 
mands are lower, the interior point methods proved su- 
perior. 



5.2 Obtaining and Evaluating Integral Solutions 

Solving the linear programming relaxation provides a 
fractional solution as well as a lower bound on the orig- 
inal integral problem. Since the fractional solutions are 
often close to integral, we derive an integral solution by 
simply thresholding the region variables. This already 
defines a segmentation, but we would also like to know 
its energy, so we have to infer the associated boundary 
variables. 

We already discussed the complications in case two 



http: //www. coin- or . org/proj ects/Clp . xml 



or more regions meet in a point (see Section 3.3). Our 
method will then select the cheapest allowed boundary 
configuration according to the selected constraint set. 
In contrast, the recent method of El-Zehiry and Grady 
[TT] will take the sum of all configurations here. 

Due to these difficulties we so far did not imple- 
ment a specialized routine to compute the energy of 
a segmentation, although this should be possible. In- 
stead, we simply re-run the linear programming solver, 
this time with all region variables fixed according to the 
segmentation. Note that this strategy was also pursued 
in [25], so the gaps reported there are w.r.t. the integer 
program, not the model itself. 

In case crossing prevention was selected we again 
add violated constraints in passes. In addition, we fix all 
"impossible" boundary variables to 0, where impossible 
refers to pairs of line segments where along one of the 
line segments the segmentation stays constant. In the 
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vast majority of cases this produced an integral solution 
and hence the optimal boundary configuration. In a few 
cases we got up to 12 fractional variables. We presently 
assume that the computed cost is close enough to the 
actual cost. 



6 Experiments 

In this section we evaluate the proposed scheme for both 
image segmentation and inpainting, where in all cases 
we consider the intensity range [0,255]. For image seg- 
mentation we also evaluate how close to the global op- 
timum we got. The experiments were run on a 3.0 GHz 
Core2 Duo machine equipped with 8 GB of memory. 
For segmentation we used the dual simplex method in 
the solver Clp, for inpainting we used the interior point 
method of Gurobi. 



Interactive Image Segmentation Next we turn to the 
problem of interactive image segmentation for color im- 
ages, where in addition to an image /:]?—>• M 3 we are 
given a set of foreground and background seed nodes as 
specified by a user. Since our focus is on evaluating the 
novel method, we did not refine the seed nodes. Hence, 
for each image the seed nodes were specified only once. 

From the given seed nodes, we estimate normalized 
histograms of color values, resulting (after smoothing) 
in distributions Pf(-) and pb(-) that are then used in 
the model: 



log(p F (/(x))[l-u(x)]dx 



log(p s (J(x))^(x) dx 



+ v\C\+\J \K C ^)\ 2 dV}(s) . 



(15) 



(16) 



6.1 Image Segmentation 

For binary image segmentation, we show experiments 
for a totally unsupervised problem and an interactive 
one where seed nodes are given. 



Unsupervised Image Segmentation For unsupervised 
image segmentation we use data functions go,g\ as in 
the piecewise constant functional of Mumford and Shah 
[19], i.e. where the intensity of an image point is com- 
pared to the mean value of a region. This results in the 
model: 

J (j(x) - /i ) 2 [l - u{x)] dx + J (j(x) - /ii)^(x) dx 

+ v\C\+\J \K C (y)\ 2 dH\y) . (14) 

c 

We set the mean values /io, \i\ to the minimal and max- 
imal intensity in the given image, respectively. 



Figure 11 shows results of our method on images 
of size 128 x 128 and 160 x 107, respectively, using an 
8-connectivity. Here we provide results for different cur- 
vature weights and it can be seen that even for very high 
weights long and thin structures are preserved. At the 
same time, the relative gap increases with the curvature 
weight. 

These results took roughly 4 hours computing time, 
where up to 9 passes were needed. Though we re-used 
the existing solution, the first passes often took as long 
as solving the initial program. 



This function is minimized over all u : ft —> [0, 1] that 
are consistent with the seed nodes. 

For the experiments, we use a 16-connectivity and 
fix the ratio of the curvature weight A over the length 
weight v to 20.0. The presented results were then ob- 
tained within half a day and using between 6 and 8 GB 
of memory. 



Figure 12 compares our results with a simple thresh- 
olding scheme and the length based method of Sec- 
tion [2] These results clearly show the tendency of 
length-based methods to suppress long and thin ob- 
jects. With curvature regularity this is remedied. 



Effects of the Constraint Sets Above we indicated that 
it is debatable whether self-intersecting region bound- 



aries should be allowed or not. In Figure [13] both pos- 
sibilities are explored. For these images there were sig- 
nificant differences, for others, like the giraffe image, 
none at all. For the boat the results clearly improve 
when forbidding self-intersections, for the other image 
they grow slightly worse. In future work we hope to get 
closer to the respective global optima to find out which 
model is actually better. 

Lastly, in this work we added the boundary consis- 
tency constraints to the original formulation in [25 . 
The latter formulation did in fact not represent the 
model correctly. Still, due to the relaxation adding the 
new constraints does not entirely solve the problem. 
When comparing the thresholded solutions of both ver- 
sions we noticed as many changes to the good as to the 
bad. In future work we hope to improve on this. 
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input 



A = 1000, v = 10 
within 0.75%. 



A = 10000, v = 10 
within 3.6%. 



A = 100000, v=!0 
within 4.4%. 




A = 100000, i/ = 10 
within 3.6%. 

Fig. 11 Evaluation of the proposed method for unsupervised image segmentation and an 8-connectivity. Shown are segmentations 
for different curvature weights and how close they are to the lower bound (and hence the global optimum). 




with crossing prevention 



without crossing prevention 



with crossing prevention 



without crossing prevention 



Fig. 13 Without crossing prevention there is a bias towards 
triple points. For the top row significantly different segmentations 
are obtained with and without the constraints. For the bottom 
row there are only minor changes. 



6.2 Inpainting 



We now turn to the problem of inpainting, where we 



use interior point solvers. Figures 14 and 15 show that 
our method is well-suited for structured inpainting, and 
that length regularity generally does not work well. 



In this work we have improved upon our work [25] 
by previously estimating the direction of incoming level 
lines and giving a tighter constraint system. Figure 16 
shows that these changes really improve the results. 



7 Conclusion 

We have presented new theory and methods for length- 
and curvature-based regularization, both for image seg- 
mentation and inpainting. For curvature (in a region- 
based context) we are the first to propose a global ap- 
proach in the sense that it is independent of initializa- 
tion. 

The results clearly demonstrate that curvature reg- 
ularity outperforms length-regularity in the presence of 
long and thin objects. Experimentally we showed that 
for segmentation our strategy of solving a linear pro- 
gramming relaxation is usually within 5% of the global 
optimum. In some cases it even finds the global opti- 
mum. 
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global optimum 




within 1.5%. 



Fig. 12 Comparison of length-based and curvature-based methods. Left column: input images with seed nodes super-imposed. Middle 
left: thresholding scheme. Middle right: length-based segmentation (proposed method). In all cases A = 1. Right: length- and curvature- 
based segmentation (proposed method) and how close the results are to the global optimum. In all cases we set A = 4 and u = 0.2. 
Images taken from the Berkeley database http://www.eecs.berkeley.edu/Research/Projects/CS/vision/grouping/segbench/. 
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